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Представлена численная методика определения дву- и трехмерного упругопластического напряженно-деформированного сос- 
тояния твердого деформируемого тела с помощью вариационно-разностного метода, реализующего вариационный принцип 
Лагранжа методом конечных разностей. Физические соотношения принимаются согласно теории малых упругопластических де- 
формаций, а геометрические соотношения берутся в виде уравнений Коши. Физически нелинейная задача решается методом 
переменных параметров упругости. На примере задач о деформировании тела вращения оживальной формы и цилиндричес- 
кого сектора проводится сравнение прямого и итерационного методов решения системы линейных алгебраических уравнений 
большого порядка, к которой приводит использование необходимого условия экстремума сеточного аналога функционала пол- 
ной потенциальной энергии системы. 


1. Введение 

Математическое описание процесса деформи- 
рования упругого тела в самом общем случае 
представляется в виде вариационных принципов 
теории упругости. Эти принципы включают в себя 


некоторые основные теоремы в форме интеграль- 
ных равенств, связывающих напряжения, дефор- 
мации и перемещения во всем объеме деформиру- 
емого тела. Вариационные принципы представля- 
ют собой теоретическую основу современных чис- 
ленных методов, позволяющих находить эффек- 
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тивное решение задач в тех случаях, когда аналити- 
ческое интегрирование уравнений теории упругос- 
ти не представляется возможным. 

Для решения упругопластической задачи вос- 
пользуемся вариационным принципом Лагранжа 
<55=0 (1) 

где Э - квадратичный функционал полной потен- 
циальной энергии системы "тело-нагрузка", кото- 
рый может быть записан для любой конструкции. 
Следствием этого является широкое использова- 
ние принципа Лагранжа при решении прикладных 
задач механики деформируемого твердого тела 
(МДТТ). В отличие от вариационных принципов 
Ху-Вашицу и Рейсснера, принцип Лагранжа явля- 
ется экстремальным, ибо для действительного нап- 
ряженно-деформированного состояния функцио- 
нал полной потенциальной энергии системы Э 
достигает минимума по сравнению с другими воз- 
можными состояниями упругого тела. Возможна 
следующая формулировка вариационного принци- 
па Лагранжа: "... из всех возможных напряженных 
состояний упругого тела, для которых выполнены 
граничные условия в перемещениях, а деформации 
и перемещения связаны геометрическими уравне- 
ниями теории упругости, истинным является то, 
при котором потенциальная энергия системы по- 
лучает минимальное значение" [1]. 

Уравнения равновесия и статические гранич- 
ные условия служат уравнениями Эйлера для 
функционала энергии Э. Поэтому при решении за- 
дачи теории упругости и пластичности в вариаци- 
онной постановке (1) нет необходимости заранее 
удовлетворять граничным условиям в напряжени- 
ях, ибо они удовлетворяются автоматически при 
нахождении минимума функционала полной по- 
тенциальной энергии системы. Вариационный 
принцип Лагранжа является средством получения 
уравнений равновесия и граничных условий в том 
случае, когда традиционный способ их выведения 
не может быть применен [2]. 

Если при минимизации функционала Э ис- 
пользуются конечно-разностные представления, 
то получаемый метод носит название вариацион- 
но-разностного метода (ВРМ). ВРМ сводит проб- 
лему минимизации функционала энергии системы 
Э, являющейся квадратичной функцией относи- 
тельно деформаций и перемещений, к задаче ми- 
нимизации функции многих переменных, отнесен- 
ных к узлам конечно-разностной сетки. 

После дискретизации функционала энергии Э 
задача об отыскании минимума функционала сво- 
дится к отысканию минимума функции многих пе- 
ременных Е (сеточного аналога функционала). 
При этом возможны два способа: прямая миними- 
зация и использование необходимого условия 
экстремума функции многих переменных, приво- 
дящего к системе линейных алгебраических урав- 
нений (СЛАУ). 


Прямые методы минимизации функции многих 
переменных позволяют находить экстремаль за- 
данной функции непосредственно, минуя уравне- 
ния Эйлера. Использование преимуществ этого 
подхода возможно при наличии эффективных ал- 
горитмов минимизации. Для прямой минимиза- 
ции функции многих переменных применяются, в 
основном, три метода: локальных вариаций, соп- 
ряженных градиентов и наискорейшего спуска. К 
недостаткам этих методов следует отнести медлен- 
ную сходимость для функций овражного типа, не- 
обходимость наличия хорошего нулевого прибли- 
жения, возможность попадания в локальный ми- 
нимум и большое время решения. 

Методы реализации минимума функции, осно- 
ванные на использовании принципа стационар- 
ности, связаны с получением и решением систем 
алгебраических уравнений высокого порядка, что 
оказывается в некоторых случаях менее выгодным. 
Эти методы можно разделить на итерационные и 
прямые. Главное отличие итерационных методов 
от прямых заключается в том, что первые дают ре- 
шение системы уравнений лишь как предел после- 
довательности итерационных приближений {/(} 
при к -» оо. Прямые методы позволяют получать 
точное решение за конечное число вычислений 
(если нет ошибок округления). Алгоритм прямого 
метода в отличие от итерационного, который сос- 
тоит из повторного применения простого алгорит- 
ма, сложен, и состоит из неповторяющихся частей. 
Метод исключения Гаусса решения СЛАУ является 
примером прямого метода. 

Представленные ниже результаты получены с 
использованием необходимого условия экстрему- 
ма полной потенциальной энергии системы. При 
этом для решения СЛАУ были опробованы оба ме- 
тода и сделаны оценки их достоинств и недостат- 
ков применительно к дву- и трехмерным задачам 
теории упругости и пластичности. 

Разработанные вариант ВРМ и численные ме- 
тодики позволяют решать статические и квазиста- 
тические упругопластические задачи в дву- и трех- 
мерной постановках на неравномерных сетках для 
тел сложной формы и сложным законом распреде- 
ления физико-механических характеристик 
(ФМХ) с целью определения их напряженно-де- 
формированного состояния и анализа прочности, а 
также создания рационального по прочности про- 
екта конструкции. 

2. Вариационно-разностный метод решения 

задачи теории упругости и пластичности 

Математическая формулировка задачи линей- 
ной теории упругости в декартовой системе коор- 
динат записывается с помощью дифференциаль- 
ных уравнений в частных производных: 

- уравнения равновесия с%+і\=0; 

- соотношения деформации-перемещения щ=(щ+ц)/ 2; 
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- соотношения напряжения-деформации сг,=%%; 

- или обратные им соотношения е^Ъ^Ои, 

- граничные условия в напряжениях <7^=7] на 8 а \ 

- граничные условия в перемещениях и,-=й) на 8 и . 

Здесь Рі - компоненты массовых сил, отнесен- 
ные к единичному объему; 7] - компоненты задан- 
ных внешних сил, отнесенные к единице площади 
поверхности; щ - заданные компоненты переме- 
щений (/ = 1,2,3); Ь і]к і - компоненты тензоров 
четвертого ранга упругих постоянных; я, - направ- 
ляющие косинусы единичной внешней нормали к 
поверхности; 8 а - часть поверхности 5” тела, на ко- 
торой заданы внешние нагрузки; 8 и - часть поверх- 
ности 8 тела, на которой заданы перемещения; 
5=5 а +5„. 

Имеем пятнадцать уравнений относительно 
пятнадцати неизвестных, а именно: шесть компо- 
нент тензора напряжений с^, шесть компонент 
тензора деформаций % три компоненты тензора 
перемещений и,-. 

Но существует другой путь решения данной за- 
дачи - вариационный. Для этого исходная задача 
теории упругости записывается в интегральной 
форме, связанной с вариационными принципами, 
и заключается в получении некоторого функцио- 
нала - интеграла от искомых функций. Таким 
функционалом является полная потенциальная 
энергия системы Э. Используя вариационный 
принцип Лагранжа (1), заменяем описанную выше 
краевую задачу вариационной задачей, в которой 
функционал энергии Э определяется классическим 
соотношением 

Э = Ѵ-А 1 -А 2 , (2) 

где II - потенциальная энергия деформации тела; 
А ь А 2 - работа объемных и поверхностных сил на 
вызванных ими перемещениях. Или, с использова- 
нием выражений для составляющих функционал 
членов 

Э = - 1 а ; . а-. йѴ - Г Р-и-сІѴ - Г Т-и-АЗ. (3) 

2ѵ У У ѵ 1 1 х 1 1 

Таким образом, решение задачи состоит в на- 
хождении поля перемещений и(х,у,і), ѵ(х,у,р), 
\ѵ(х,у,%) в направлении осей координат х, у, і соот- 
ветственно, доставляющих минимум функционалу 
полной потенциальной энергии системы Э. 

Одним из способов реализации минимума 
функционала полной потенциальной энергии сис- 
темы является вариационно-разностный метод. 
Сущность ВРМ, реализующего вариационный 
принцип Лагранжа (1) с помощью метода конеч- 
ных разностей, заключается в сведении задачи ми- 
нимизации функционала полной потенциальной 
энергии Э, являющейся квадратичной функцией 
относительно деформаций или перемещений, к за- 
даче минимизации функции многих переменных, 
отнесенных к узлам конечно-разностной сетки. 
Сеточный аналог Е(и„ѵ р ,щ) функционала (3) реа- 


лизуется заменой составляющих его интегралов ко- 
нечными суммами, а производных - их конечно- 
разностными представлениями. Приближенное 
конечно-разностное выражение полной потенци- 
альной энергии системы Е получается в предполо- 
жении, что в пределах каждой ячейки все функции 
и их производные остаются постоянными. Следует 
подчеркнуть достоинства вариационно-разностно- 
го метода, основанного на вариационном принци- 
пе Лагранжа: простота математической формули- 
ровки задачи; ясный физический смысл использу- 
емого функционала; автоматическое выполнение 
уравнений равновесия и статических граничных 
условий; возможность использования метода для 
расчета тел сложной формы, в том числе неодно- 
родных по деформационно-прочностным характе- 
ристикам материалов; сводимость проблемы к ре- 
шению системы линейных алгебраических уравне- 
ний, для реализации которой существует достаточ- 
но надежный математический аппарат линейной 
алгебры. Причем матрица системы имеет ленточ- 
ную структуру, симметрична и положительно опре- 
делена, что очень важно при численной реализа- 
ции для практических приложений; возможность 
поэтапного контроля точности выкладок и проме- 
жуточных результатов. 

Одной из самых первых работ, где был приме- 
нен вариационный подход для получения сеточных 
уравнений при численном расчете изгиба прямоу- 
гольных армированных пластин поперечной наг- 
рузкой, является работа Л. А. Оганесяна [3]. Анало- 
гичный подход к решению задачи о деформирова- 
нии пластин рассмотрен в [4, 5]. Применительно к 
решению осесимметричных и плоских задач тео- 
рии упругости ВРМ был разработан Д.С. Гриффи- 
ном и Р.Б. Келлогом [6]. В отличие от выше упомя- 
нутых работ, где реализация СЛАУ проводится ите- 
рационным методом и методом сопряженных гра- 
диентов, здесь используется метод Гаусса. 

В общем случае функционал энергии является 
неквадратичным. Для анализа упругопластическо- 
го поведения конструкций используются либо де- 
формационная теория пластичности А.А. Илью- 
шина, которая хорошо описывает процесс упру- 
гопластического деформирования при монотон- 
ном нагружении и широко используется в расчетах 
упругопластических тел, или теория идеально- 
пластического тела (теория течения), устанавлива- 
ющая связь между приращениями напряжений, де- 
формаций и некоторыми параметрами пластичес- 
кого состояния. Ниже используется деформацион- 
ная теория пластичности. Решение физически не- 
линейной задачи сводится к последовательности 
решений квадратичных задач с уточняемыми в 
каждом приближении параметрами нелинейности. 
С помощью этих параметров и осуществляется 
движение вдоль диаграммы нелинейной зависи- 
мости интенсивности напряжений от интенсив- 
ности деформаций с г, ~ е„ которая аппроксимиро- 
валась ломаной двухзвенной линией. 


25 




Известия Томского политехнического университета. 2003. Т. 306. № 3 


Изложенный алгоритм неоднократно был про- 
тестирован на многочисленных примерах решения 
упругих и упругопластических задач, а также срав- 
нением с экспериментом. С помощью созданного 
варианта ВРМ был решен ряд статических и ква- 
зистатических упругопластических задач для суще- 
ственно неоднородных тел вращения в дву- и трех- 
мерной постановках. 

Используя интегральную формулу преобразова- 
ния объемного интеграла в поверхностный интег- 
рал, запишем выражения для производной по объ- 
ему от скалярной функции Г в декартовых коорди- 
натах х, у, і [7]: 

дР / дх = Ііш |[> Рп х сІ5 / V, 

дР / ду = Ііш {_[ РПуСІ8 / V, 

дР Іді = \\т§Рп 2 с18 IV, 

где п х , п у , п- - косинусы углов между вектором нор- 
мали п к замкнутой поверхности 8, заключающей в 
себя объем V, и соответствующей осью координат: 

п х - С08(п, І), П = С08(п, І), П. = С08(п, к). 

Таким образом, пространственные производ- 
ные в трехмерном пространстве определяются че- 
рез интеграл по замкнутой поверхности. При этом 
функция С должна быть непрерывной в объеме V и 
иметь там непрерывные ограниченные первые 
частные производные. Накладываются некоторые 
ограничения и на поверхность 5. 

Аналогичные формулы имеют место для дву- 
мерного осесимметричного случая (так называе- 
мые ноховские производные) [8]: 

дР / дг = Ііш } Рп г с18 / А, 

Я->0 С 

дР / дг = Ііш [ Рп 2 сІЗ / А, 

А->0 С 

где С - граница области А\ 8 - длина дуги; п„ п. - 
косинусы углов между вектором нормали п к внеш- 
нему контуру С и осями цилиндрической системы 
координат г, ср, і. 



Рис. 1 . Шаблон для случая дискретизации расчетной облас- 
ти с помощью четырехугольных ячеек 


Дискретизация расчетной области в двумерной 
задаче проводится с помощью четырехугольных 
ячеек (рис. 1), а в пространственной трехмерной 
задаче - с помощью многогранников (шестигран- 
ников), ограниченных четырехугольными плос- 
костями (рис. 2). В частном случае это параллеле- 
пипеды (прямоугольные в том числе) и призмы с 
основаниями в виде трапеций и с параллелограм- 
мами и прямоугольниками в качестве граней. 



Рис. 2. Шаблон для случая дискретизации расчетной облас- 
ти с помощью многогранников 

После дискретизации вариационной задачи и 
получения выражения сеточного аналога функцио- 
нала энергии использование необходимого усло- 
вия экстремума 

дЕ(ц„ѵ р ,уѵ к )/ди, =0, дЕ(и, ,ѵ р ,щ )/ дѵ р =0, 
д Ш(и 5 ,ѵ р ,\ѵ к )/д\ѵ к =0, 

(5 = ^,5! +1,..., 5; р = р 1 ,р 1 + 

к = + 1 , . . . , К\ 

сводит задачу минимизации функции многих пере- 
менных Е(и„,ѵ р ,щ) к решению системы п линейных 
алгебраических уравнений относительного п иско- 
мых компонент вектора перемещений узлов конеч- 
но-разностной сетки 

ХЧ/у =*/>* = !» 2, ... ,п. (4) 

}= і 

Здесь коэффициенты щ - элементы матрицы 
{А} СЛАУ; Д - проекции и„ѵ р ,щ вектора перемеще- 
ний на оси координат; Ъ, - свободные члены, вклю- 
чающие в себя статические и геометрические гра- 
ничные условия. Следует отметить, что в задачах 
МДТТ матрица коэффициентов {А} системы урав- 
нений (4), к которой приводит алгоритм ВРМ, 
симметрична, положительно определена и имеет 
ленточную структуру [4]. 
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Естественные науки 


При получении аналитического выражения ус- 
ловия стационарности функции энергии Е для оп- 
ределения коэффициентов а ѵ необходимо пользо- 
ваться шаблонами, представленными на рис. 1,2. 
Структура выражения функции энергии такова, 
что от смещения в некотором узле сетки с номером 
О зависит лишь часть энергии Е, а именно энергия 
ячеек, примыкающих к этому узлу. Количество та- 
ких ячеек колеблется от одной (угловой узел) до че- 
тырех и восьми ячеек (в середине массива) для дву- 
и трехмерной задачи соответственно. 


3. Итерационные методы решения СЛАУ 


Из итерационных методов были рассмотрены 
методы Зейделя (Гаусса-Зейделя) и верхней релак- 
сации (МБР), которые получили наибольшее расп- 
ространение при итерационном способе реализа- 
ции СЛАУ, а также МБР с оптимальным коэффи- 
циентом релаксации (ОКР) <о т . Алгоритм метода 
Зейделя, сходимость которого доказана на основе 
общей теории разностных схем [9], в предположе- 
нии, что диагональные элементы матрицы отлич- 
ны от нуля (а„^0), записывается в следующем виде: 


7=1 


(*+ И 


+ .?У/ )= і. 

у=г+ 1 


где //'-./-я компонента вектора решения итераци- 
онного приближения с номером к. 

Величина перемещения ^на (&+1)-ой итерации 
определяется из уравнения 

7=1 7=1+1 

Из последнего соотношения видно, что метод 
Зейделя относится к так называемым неявным 
двухслойным итерационным методам. 

Для ускорения сходимости метода Зейделя его 
модифицируют введением в итерационную схему 
коэффициента релаксации со и приходят к так на- 
зываемому методу верхней релаксации. Этот метод 
был предложен Франкелом [10] и Янгом [11] и 
обобщен другими авторами. Основные результаты 
изложены в работе Варги [12]. В этом случае итера- 
ционный процесс осуществляется по следующим 
формулам: 

7=1 

-ю Е а„/} к) +а)Ь г (5) 

7=1+1 

При этом реализация одного итерационного 
шага осуществляется примерно тем же количест- 
вом действий, что и в методе Зейделя. Этот моди- 
фицированный метод при со > 1 носит название ме- 
тода верхней релаксации, при со = 1 - полной ре- 
лаксации (метод Зейделя) и при со < 1 - нижней ре- 
лаксации. Доказано, что для 0 < со < 2 метод всегда 


сходится и расходится при со < 0 или со > 2 (при 
со = 0 или со = 2 не сходится и не расходится) [12]. 
При этом матрица коэффициентов {А} должна 
быть симметричной, неособенной, положительно 
определенной, а все диагональные коэффициенты 
СЛАУ а,і > 0. Оптимальное значение со 0р , коэффи- 
циента со, обеспечивающее наиболее быструю схо- 
димость итерационного процесса, лежит в преде- 
лах 0 < со < 2. 


Рассмотрим вопрос об определении ОКР в ме- 
тоде верхней релаксации, когда 1 < со < 2. В работе 
[13] приведена следующая формула для определе- 
ния С0 ор ,\ 



( 6 ) 


где Р =Вт 


у(А+1) _у( 


<‘->11 




. 11/КЕё 2 , к-\,к, к+1- 


номера итераций, выполняемых при со = 1. 

На практике для вычисления р 2 и со ор , делается 
достаточно большое число итераций при значении 
со = 1 , а затем решение продолжается при со = со 0рІ , 
определенным из (6). 

А.П. Деругой и А.А. Ларионовым [14] для уско- 
рения сходимости итерационного процесса (5) 
предложена модификация способа (6) определения 
со 0р , в виде следующей зависимости от числа итера- 
ций к: 

со орІ (к) = а + Ье- ск . (7) 


Чтобы определить коэффициенты а, Ъ, с, для за- 
данных величин к = / - А, I, I + А вычисляются по 
формуле (6) значения со„ р , и приравниваются выра- 
жению (7): 

со орІ (/ - А) = а + Ъе~ с °~ А) , со ор1 (/) = а + Ъе с1 , 
(о ор1 (1 + А) = а + Ъе- с(1+А> . 

Из решения полученной системы уравнений 
получаем 

_ п\ і Р\Рі 

а ®орі СО 3" •> 

р Р\ — Р 2 


где Р\ (Оіірііі) сО/ )р! ( I А) ■ Рі со о А I + А) со //р! ( I) . 

Для к = о о формула (7) для ОКР запишется в 
следующем виде: 

Рі Рі 

Окончательно, ввиду того, что Ъет со ° = 0, имеем 

а >о Р Л с °) = °>о Р Л 1 ) + ~Ріуг- (8) 

Р\ Рі 

Таким образом, вычислив ОКР по формуле (8), 
итерационный процесс (5) продолжаем со значе- 
нием аѵ(°°). 
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